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Hydration of hydrophobic solutes in water is the cause of different phenomena, including the 
hydrophobic heat-capacity anomaly, which are not yet fully understood. Because of its topicality, 
there has recently been growing interest in the mechanism of hydrophobic aggregation, and in the 
physics on which it is based. In this study we use a simple yet powerful mixture model for water, 
an adapted two-state Muller-Lee-Graziano model, to describe the energy levels of water molecules 
as a function of their proximity to non-polar solute molecules. The model is shown to provide an 
appropriate description of many-body interactions between the hydrophobic solute particles. The 
solubility and aggregation of hydrophobic substances is studied by evaluating detailed Monte Carlo 
simulations in the vicinity of the first-order aggregation phase transition. A closed-loop coexistence 
curve is found, which is consistent with a mean-field calculation carried out for the same system. 
In addition, the destabilizing effect of a chaotropic substance in the solution is studied by suitable 
modification of the MLG model. These findings suggest that a simple model for the hydrophobic 
interaction may contain the primary physical processes involved in hydrophobic aggregation and in 
the chaotropic effect. 



I. INTRODUCTION 

Many recent studies have focused on the surprising 
thermodynamic properties associated with the hydration 
of hydrophobic substances in aqueous solutions 0, @- 
Indeed, the heat capacity of non-polar solutes in water 
is large and positive at room temperature. This behav- 
ior, known as the "hydrophobic heat-capacity anomaly" , 
stands in sharp contrast to the observations made for hy- 
drophilic solutes in water Q. Experiments confirm that 
the solubility of small hydrocarbons in water decreases 
when the solution is heated near room temperature. It 
is generally believed that these thermodynamic phenom- 
ena are associated with a structural change in the sol- 
vent, and the hydrophobic effect is therefore considered 
as directly related to the anomalous properties of liquid 
water p|. 

Due to hydrogen bonds between molecules, pure water 
is highly structured. Adding a non-polar solute, which 
is unable to form hydrogen bonds, disrupts the struc- 
tural arrangement by breaking bonds, thus raising the 
enthalpy of the solution and increasing considerably the 
entropy. Because the entropy change dominates over the 
enthalpy change at room temperature, the hydration pro- 
cess is considered to be "entropy driven" . 

Although this intuitive interpretation provides a first 
understanding of the hydrophobic effect, it still suffers 
from a number of limitations. In contrast to the expected 
increase in enthalpy due to the breaking of hydrogen 
bonds, the process of transferring a non-polar particle 
to water from the liquid state in pure solute is found ex- 
perimentally to be enthalpically even slightly favorable 
(negative), at least at low temperatures This phe- 
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nomenon can be explained only by a rearrangement of 
the water molecules in the vicinity of the solute, result- 
ing in a recovery of the lost hydrogen bonds which tend to 
be slightly stronger than before This new structure, 
however, causes an entropy loss attributed to an increase 
in local order. Measurements of the entropy and enthalpy 
changes for different solute sizes suggest that this local 
ordering of water molecules around the hydrophobic par- 
ticle is not unique, but that a number of different orga- 
nizations is possible 0. Hence, the anomalous thermo- 
dynamic properties are a consequence of changes in the 
state of water molecules due to the insertion of hydropho- 
bic solute particles, rather than being accounted for by 
water-solute interactions. 

At room temperature the solubility of small, hydropho- 
bic solute species in water increases when the system 
is cooled. At sufficiently low temperatures a homoge- 
neous mixture is found, provided that none of the com- 
ponents crystallizes. Heating the mixture produces a de- 
crease in solubility, which can result in a phase transition 
to a state with two phases of different solute densities. 
This transition temperature is called the Lower Criti- 
cal Solution Temperature (LCST), and has been mea- 
sured for different solutions of hydrophobic particles in 
water 0,0,0,^3- After reaching a minimum, the solubil- 
ity increases steeply at higher temperatures. If a coexis- 
tence of two phases is found near room temperature, an- 
other phase transition to one homogeneous phase occurs 
at the Upper Critical Solution Temperature (UCST), on 
condition that the boiling point of the solution is not 
reached [Tlj . 

Aqueous solutions of non-polar particles which show 
a LCST also have an UCST and therefore a closed-loop 
coexistence curve in the phase diagram, provided that 
none of the liquids undergoes a phase transition to a 
gaseous phase before that temperature is reached. This 
occurs because the rising temperature increases the en- 
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tropy, which favors a homogeneous mixing of the com- 
ponents. In addition, the number of hydrogen bonds 
formed, which is responsible for the phase separation, de- 
creases with increasing temperature. This type of closed- 
loop miscibility curve has been found in different binary 
solutions including nicotine/water and poly (ethylene gly- 
col) /water 0,0; as well as by modeling of binary mix- 
tures containing hydrogen-bonded molecules |l2|, [l3|, [l4| . 

A complete understanding of the molecular mecha- 
nisms underlying the hydrophobic effect is essential in 
order to explain a variety of phenomena, including the ag- 
gregation of hydrophobic solutes in water and the desta- 
bilization, denaturation and aggregation of proteins, pro- 
cesses which are responsible for many life-sustaining pro- 
cesses but also for different diseases 0, E3, • 

Aggregation of non-polar particles in aqueous solu- 
tions can be attributed to the effective hydrophobic in- 
teraction between solute particles. This interaction is 
mainly solvent-induced, because the approach of two hy- 
drophobic solute particles reduces the surface exposed 
to water compared with the total surface when they are 
completely separated. However, despite its importance 
for the understanding of a number of still-unexplained 
phenomena, the study of hydrophobic aggregation is at 
present incomplete and needs further investigation. In 
this work we describe and analyze solutions of hydropho- 
bic particles, particularly in the vicinity of the aggrega- 
tion phase transition, using a simple two-state model in 
three dimensions. The results of a mean-field calculation 
are compared to Monte Carlo simulations. In addition, 
the model is extended to include the effect of a chaotropic 
substance in the solution, which causes an increase in sol- 
ubility of the non-polar solute for reasons which are at 
present not well understood. 

This paper is organized as follows. Section II intro- 
duces the model used for a detailed description of the 
solution. In Section III a mean-field calculation is per- 
formed, giving a first idea of the phase diagram near the 
aggregation transition. Monte Carlo simulations confirm- 
ing these findings are presented in Section IV. In Section 
V the model is adapted to describe the chaotropic effect. 
Finally, Section VI offers a discussion of the results, as 
well as some suggestions for further investigations. 



II. MODEL 

In order to describe the interaction between hydropho- 
bic solutes and water, a simple model which describes the 
essential physics of the system is required. It is generally 
believed that the driving force in the aggregation pro- 
cess is the effective repulsive hydrophobic interaction be- 
tween the polar water and non-polar solute . As early 
as 1945, Frank and Evans noted that this interaction 
arises from a cage-like arrangement of water molecules 
around the solute, which allows them to optimize their 
mutual hydrogen bonding and thus to minimize their en- 
ergy with respect to bulk liquid water. Entropically, how- 



ever, this strict ordering is unfavorable compared with 
the disordered configurations in bulk water. 

The minimal model which contains the essential fea- 
tures of the physics of water as an aqueous solvent is 
the bimodal description of Muller, Lee, and Graziano 
(MLG) [SI13. Thc model describes liquid water by di- 
viding it into two different populations based on the num- 
ber of hydrogen bonds formed. Water molecules which 
are highly hydrogen-bonded to their neighbors have fewer 
rotational degrees of freedom and thus a lower multi- 
plicity of degenerate configurations (lower entropy) than 
unbonded molecules. These are therefore denoted as "or- 
dered" , while water molecules with many broken hydro- 
gen bonds are considered to be "disordered". The pres- 
ence of a non-polar solute alters the enthalpy and entropy 
of water molecules in the solvation shell of each solute 
particle, and so a further distinction is required between 
"shell" and "bulk" water. More details of the physical 
processes underlying these considerations are provided 
below. This simple model has been justified recently from 
molecular simulations of water [2l| . 

In this form the MLG model reproduces correctly the 
ordering and strengthening of the hydrogen bonds in the 
first solvation shell of an added non-polar solute molecule 
at low temperatures, as well as the opposite behavior at 
high temperatures. It has been shown that the model 
provides an adequate description of the heat-capacity 
anomaly [2^, and it has been used to reproduce con- 
sistently the important properties of protein solutions, 
including warm and cold denaturation |23ll24j. 

In this analysis we use an adapted version of the MLG 
model. An appropriate description of the solvent in the 
vicinty of hydrophobic solute particles which may be 
much larger than individual water molecules is obtained 
by allowing each site of a lattice representing the sys- 
tem to be occupied by a group of water molecules. The 
bimodal nature of the MLG model is preserved in this 
coarse-grained version by specifying only two types of 
water cluster at each site, where now an "ordered" site 
is characterized by having most of the hydrogen bonds 
among the molecules in the cluster intact, while a "dis- 
ordered" site is understood to have a number (but by 
no means all) of these bonds broken. While this is the 
simplest possible approximation to thc continuous dis- 
tribution of intact or broken bonds within a site, we will 
demonstrate that it retains the capability of describing all 
the primary physical properties observed in aqueous so- 
lutions of non-polar solutes. The energy and degeneracy 
parameters for the coarse-grained model are determined 
by the same processes as in the bare MLG description 
outlined above, which we now discuss in more detail. 

An approximation where a group of water molecules is 
considered as one entity is justified when the non-polar 
solute particle is relatively large compared with a single 
water molecule. In this case the formation of a complete 
cage around the solute particle is rather improbable be- 
cause it must be formed rapidly in the presence of local 
thermal fluctuations, and may even be prevented steri- 
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cally |25j. Partial cages may therefore be formed in the 
vicinity of a solute particle, rather than one complete 
cage. In addition, formation of a hydrogen bond pro- 
motes the formation of further hydrogen bonds, which 
are stronger than before due to the change in charge dis- 
tribution on forming the first bond 26]. This mutual 
reinforcement is known as "cooperativity" of hydrogen 
bonds, and leads to the formation of chains or clusters of 
hydrogen-bonded water molecules, whose extent depends 
on the size and the shape of the solute particles. 

Water molecules which participate in hydrogen-bonded 
clusters have a higher degree of order and fewer rota- 
tional degrees of freedom than those in regions with 
many unbonded molecules 25] . The number of possi- 
ble configurations of such an "ordered" cluster is thus 
significantly smaller than that of a "disordered" group 
of water molecules for both shell and bulk water. For 
steric reasons, fewer hydrogen-bonded water configura- 
tions arc possible around a non-polar solute particle 
which is unable to form hydrogen bonds. These shell 
water molecules are forced into a tangential orienta- 
tion |2jj, whereas the molecules in bulk water may also 
form radially oriented hydrogen bonds with central wa- 
ter molecules replacing the solute |2^. The degeneracy, 
or total number of configurations, of a hydrogen-bonded 
cluster of shell water ("ordered shell") is consequently 
smaller than that of a hydrogen-bonded cluster of bulk 
water ( "ordered bulk" ) . 

In contrast, fewer orientational configurations exist for 
unbonded water molecules in the bulk than next to a non- 
polar solute particle . The geometrical reason for this 
result is that for a shell site no hydrogen bonds are pos- 
sible in direction of the solute particle, unlike the bulk 
situation obtained on replacing the solute particle by wa- 
ter. All orientations in which water molecules form radial 
hydrogen bonds with the central water in the bulk case 
(contributing to the ordered bulk degeneracy) are there- 
fore transformed into configurations with many broken 
hydrogen bonds when the central water is replaced by a 
non-polar solute particle. We take such sites to be "dis- 
ordered" in the bimodal sense discussed above. The de- 
generacy of a group of water molecules with many broken 
hydrogen bonds is then higher in the shell ("disordered 
shell" ) than in the bulk ( "disordered bulk" ) . 

In summary, these considerations lead to a distribu- 
tion of the total number of states q according to the 
sequence qd s > qdb > qob > qos- Because many fewer 
configurations have intact than broken hydrogen bonds, 
the difference in degeneracy between ordered and dis- 
ordered states is much higher than that between shell 
and bulk states for both types of site. In fact the dif- 
ference in degeneracy between shell and bulk states de- 
pends primarily on the number of possible radial hydro- 
gen bonds, which is much smaller than the total number. 
In much of what follows we employ the degeneracy fac- 
tors q ds = 49, q d b = 40, q ob = 10, and q os = 1. These 
relative values are chosen to be qualitatively represen- 
tative of the degeneracies expected for a system of small 



solute particles in water, with both experimental and the- 
oretical [52] justification based on the above properties 
of water. We note here that, while it is known that a 
number of ordered states exists even for the cage-like wa- 
ter structures coordinating dissolved solute particles, the 
absolute multiplicities of the q factors contribute only an 
additive constant to the free energy and are irrelevant 
for the phenomena to be discussed below; thus we set 
the lowest degeneracy to q os = 1. 
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FIG. 1: Energy levels in the MLG model for water. 

Fig. shows the schematic energy levels of a group of 
water molecules in the coarse-grained MLG-type model. 
The hydrogen-bond energy is optimized in the ordered, 
cage-like shell structure where strong, tangentially ori- 
ented hydrogen bonds are formed among water molecules 
in the cluster ("ordered shell"), and hence E os is lowest. 
Direct experimental evidence for the tangential orienta- 
tion of hydrogen bonds between water molecules in the 
solvation shell of non-polar solute particles has been pro- 
vided recently in Refs. [27l |28| . When the solute par- 
ticle is replaced by water, radial hydrogen bonds may 
be formed. However, clusters including radially oriented 
water molecules have on average a higher hydrogen-bond 
energy than those whose bonds are predominantly tan- 
gentially oriented, because for steric reasons a radial hy- 
drogen bond precludes another good tangential hydrogen 
bond to a first-shell neighbor. This result was demon- 
strated in a model including oriented hydrogen bonds 
between water molecules |22(. In bulk water, both con- 
figurations are possible, and thus the average energy in 
a cluster of hydrogen-bonded bulk water E {, ("ordered 
bulk") is higher than E os . 

The energies Ed s ( "disordered shell" ) and Edb ( "disor- 
dered bulk" ) are relatively much higher than those of the 
respective ordered states, because breaking of hydrogen 
bonds is required in forming the disordered states. The 
average hydrogen-bond energy of a group of disordered 
shell water molecules decreases when the solute particle 
is replaced by water, because some radial hydrogen bonds 
broken by the solute may be formed, lowering the average 
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energy of the group and ensuring that Edb < Eds- 

Specific determination of the energy values for a se- 
lected binary system would require structural calcula- 
tions and molecular dynamics simulations p^ |. However, 
such refinement is not necessary for the general phenom- 
ena to be illustrated below. We note here that the tem- 
perature scale /3 _1 is defined by the energy scale. In the 
following we use the parameter values Eds = 1.8, Edb = 
1.0, E b = —1.0, and E os = —2.0, which are thought to 
be quite generic for aqueous systems, and furthermore 
agree closely with the energies used in a successful de- 
scription of the thermodynamic behavior of biopolymers 
in water |2^. While the results of the calculations to 
follow are not particularly sensitive to the exact values 
selected, it remains possible to refine these parameters by 
comparison with experiment to obtain good qualitative 
agreement with measured quantities for different solu- 
tions. 

On a cubic lattice the energy of a system of N sites, 
occupied either by particles (rii = 0) or by water (m = 1), 
is given by the Potts-like Hamiltonian 

N 

1=1 

+ (E bSi,a ob + Edb5i,cr db )\i], (1) 

where Ai is I if site i is surrounded only by water, is 
otherwise, and is defined as the product of the nearest- 
neighbor factors, Ai = Y[/j j\ n j- Each water site i can be 
in one of the q different states which are divided among 
the 4 energy levels shown in Fig. ^ Therefore, S^ aos is 
1 if site i is occupied by water in one of the q os ordered 
shell states and otherwise, and Si,cr da is 1 if it is occupied 
by water in one of the qds disordered shell states and 
otherwise. The same holds for the bulk states. 

An important observation is that this model does not 
include interactions between different water sites. Hence 
the model is valid as long as water is liquid, and ne- 
glects long-range effects arising from extended hydrogen- 
bonded networks. 

Hamiltonian leads to the partition function 

Z= e-? H ^>^\ (2) 

where every term of the sum represents the statistical 
weight of the corresponding configuration = fcgT) . 
Performing the sum over the configurations of the states 
of each water site, including the number of states of the 
respective energy levels, gives 

Z = ^JJ e -/3[B» i A 4 +Sn 4 (l-A 4 )] ) (3) 

where we define 

S = ~H<lo S e-e E °« + q ds e-0 E «>] and (4) 

B = -^-\n{q ob e-P E °> + q db e~^\. (5) 



The formal method for obtaining an effective Hamilto- 
nian (and effective interactions, see below) is to integrate 
over the degrees of freedom of the particles thought to 
be responsible for the interactions, which are the solvent 
molecules. The canonical partition function, i.e. the par- 
tition function for a fixed number of particles N p (whence 
the number of water sites is N w = N — N p ), is 

Z Na = J2e- 0H ^ n ^, (6) 
{«,} 

where the effective Hamiltonian -ff c ft is formally the 
free energy of the solvent at fixed particle configuration. 
In a system where the number of particles is not fixed, 
a chemical potential /i associated with the replacement 
of water sites by solute particles (/i is thus the differ- 
ence between the chemical potentials of water and solute 
particles) is included, and the grand canonical partition 
function may be expressed as 

3 = "£^Z NiB (7) 
= £y/»fl£K»«H W here (8) 

N 

H !s[{^}} = ^[(S + MK + CS-SKAi]. (9) 

i=l 

i/^[{rii}] represents an effective Hamiltonian for single 
sites, and provides the first step in obtaining the effec- 
tive interactions between particles. In order to describe 
the consequences for hydrophobic particles in the solu- 
tion, we replace the number of water molecules rii by the 
number of particles m.,, 

rii = 1 - mi, (10) 

where m-i is 1 if the site is occupied by a particle and 
otherwise. With this substitution, Ai becomes the prod- 
uct over the nearest neighbors, Ai = i)(l — 171 j)> an d 
takes the value 1 only if site i is completely surrounded 
by water molecules, or otherwise. Introducing an effec- 
tive interaction J G ff = (£> — S) between particles, and the 
effective chemical potential /i e ff = (S + /z), the effective 
Hamiltonian for the particles becomes 

N 

H !s, P [{™J] =K- ^[m™. + Jeff (mi - 1)A,], (11) 
»=i 

where K = N(S + fx). In this formulation it can be 
seen that the interactions are not limited to two-body 
terms, but include many-body interactions through the 
last term. 

Because S is negative and decreases continuously, fi c s 
decreases as temperature increases, and the solubility 
drops. If fi > —S, /i c ff is positive at low temperatures 
and the solubility is high. In the ground state (T = 0), 
this results in a solid phase for /i > —S and pure water 
for [i large and negative. 
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At low temperatures, J c g is positive and therefore the 
interaction term is minimized by = 1, which means 
that there is a repulsive force between particles. At high 
temperatures, however, increasing entropy effects cause 
J e ff to become negative, and the minimal interaction en- 
ergy is obtained for A» = 0, resulting in an attractive 
force between particles. These different effective forces 
result from the interplay of entropic and enthalpic ef- 
fects, and give rise to the complex properties of solutions 
of hydrophobic particles. 



III. MEAN-FIELD CALCULATION 

We model the system on a three-dimensional square 
lattice, where every site has z = 6 nearest neighbors, 
which can be occupied either by a particle or by water. 
Our first approach is based on the assumption that spa- 
tial fluctuations of the density are insignificant. In order 
to make qualitative predictions, a mean-field approxima- 
tion is used for the average occupancy of a site, (rii) = p, 
thus treating the density as constant throughout the sys- 
tem. (Because of the nature of the model, p is connected 
to the number of water molecules, and therefore the so- 
lute particle density, which is used in the graphs to follow, 
becomes p p = I — p.) The grand canonical mean-field free 
energy of a particle is given by 

z+l 



/ = (B - S)p z+1 + (S + p) P 
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FIG. 2: Differences in free energy F, enthalpy U, and entropy 
S (multiplied by T) per particle between a system of particle 
density p p = 0.1 and pure water. The system with density p is 
forced into a homogeneous state in which the particles are in 
solution. Between the critical temperatures Tl and Tu, pure 
water is energetically favorable and the system separates into 
two phases. 

Comparison of the free energy per particle of a sys- 
tem with a given density to that of pure water (with the 
same number of water sites) illustrates the thermody- 
namic properties of the system. Because the density is 



fixed (p = 0), the system is forced into the homogeneous 
state although the free energy of a mixture of two sepa- 
rated phases with different densities might be lower. Pure 
water with the same number of water sites has approxi- 
mately the same total separation of water and particles 
in a system of the same overall density. As shown in 
Fig. |21 a homogeneous mixture of density p p — 0.1 is en- 
ergetically favorable compared to pure water at low and 
at high temperature, indicating good solubility of the so- 
lute in water. In between, however, the solute prefers not 
to dissolve, and two phase transitions occur at the tem- 
peratures Tl and Tjj where the differences in enthalpy 
and in entropy (multiplied by T) are equal. 

Because the system is forced into a homogeneous state, 
and p — in the calculation at constant p p , the resulting 
free energy is larger than for a system minimizing f(p) at 
a finite value of p chosen to produce the same density p p . 
In addition, the comparison with pure water neglects the 
small water-solute surface of a system in which the two 
species are completely separated, as it takes into account 
only the free energy of the water phase. The immiscible 
region therefore increases for a system allowed to sepa- 
rate into two phases. All of the following calculations 
are carried out by minimizing f(p) in Eq. I|12[l. which 
provides the best available approximation. 




FIG. 3: p-T phase diagram for an aqueous solution of non- 
polar solute particles obtained by mean-field calculation. The 
outer line represents the closed-loop coexistence curve, while 
the inner line marks the spinodal curve. The arrows of de- 
creasing T demarcate quenches into the metastable (Ti) and 
spinodal (T2) regions (see text). 

Based on these equilibrium values, a phase diagram for 
the values of p can be obtained as a function of tempera- 
ture, and is displayed in Fig. [21 The outer line represents 
the closed- loop coexistence curve T co (p), outside which 
the system is in a homogeneous state. Because we find 
aggregation by heating at low temperatures, and we ex- 
pect that at high temperatures the entropy of solvation 
should dominate, a closed-loop solubility curve showing 
a LCST and an UCST [H HJ H3 is to be anticipated. 
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The inner line is the spinodal curve T sp (p). 

Microscopically, the process of phase separation may 
be observed when quenching the system from the homo- 
geneous phase at To > Tjj, at a constant density po; hito 
the coexistence region (Ti or T2). At long times after the 
quench, the system will be completely separated into two 
phases (with densities p c \ and p c ^) in a ratio depending 
on the initial density p$. The fraction of the total volume 
occupied by phase i is Vi = ^Z 9 ^. > where j stands for the 
other phase. 




FIG. 4: Comparison of p-T phase diagrams for hydrophobic 
particles in water obtained by mean-field calculation (solid 
line) and by Monte Carlo simulations (o). Both coexistence 
curves agree within the simulation error below the UCST of 
the Monte Carlo simulation T{j MC . The mean-field calcula- 
tion results in a higher UCST (Ty mf ) than do Monte Carlo 
simulations. A crystal phase appears at low temperatures in 
Monte Carlo simulations. The lines show a system heated 

at constant particle density p v = l/(z + 1) ( ) and at 

p p — 0.5 (•••)! corresponding to the vertical lines in Fig. |S] 
Although the critical density, pZ, = l/(z + 1), determined by 
mean-field calculation differs slightly from the value found by 
Monte Carlo simulations, it lies within the error, and no vis- 
ible deviation is found from the expected behavior. 

Fig. 0] shows the phase diagram of the system a func- 
tion of p and T. On heating at constant p (verti- 
cal lines) a phase transition is found at a temperature 
T£ < T t (p) < (for p* L < p < p^), where the par- 
ticle density of the system jumps discontinuously from 
a value p c i to p c \. At a constant temperature, such as 
T\ in Fig. [31 the free-energy density f(p) is minimal at 
Pd and p C 2- Heating the system from below the transi- 
tion temperature at fixed chemical potential p results in 
a discontinuous jump in density from p c \ to p C 2 at T t . 

The mean-field p-T phase diagram in Fig. 0] shows a 
first-order transition line bounded by two critical points 
characterized by the same critical solvent density p* . An- 
alytically, p* may be obtained by imposing that the two 
local minima of the free energy /, as well as the inflection 
points of /, coincide at the critical points p*, T* and p* . 
The first, second and third derivatives of the free energy 
with respect to the density must therefore vanish at the 



critical points. We calculate simultaneously 
(z + l)(B-S)p z + (S + p) 



df 

dp 



+(3- 1 Ht JL ~ ) = 0. 



1 



P 



d 2 p 



d 3 P 



= z{z + l){B-S)p z 



+/T 1 - 



f 



= 0, and 



P(l - P) 
z{z + l){z-l){B-S)p z - 2 



(i-pY p?' 



0, 



and simplify Eq. (|14|l to the form 

/r 1 = -z{z + \)(\-p)(B-S) P z 



(13) 



(14) 



(15) 



(16) 



Introducing Eq. (|16|) into Eq. (|15fl leads to the critical 
density p* — ^-j- (i.e. critical particle density p* = jxj), 
which depends only on the effective coordination number 
z. Inserting p* in Eq. HHJ) provides the LCST T£ = 
0.518 and the UCST T£ = 2.79 for the parameter values 
chosen. Finally, Eq. i|13|) gives the corresponding lower 
and upper critical chemical potentials p* L = 1.69 and 
plj = 7.53. These values are shown in the p-T phase 
diagram in Fig. 0] 



IV. MONTE CARLO SIMULATION 

To analyze the behavior of the system in more detail, 
we study a three-dimensional system of N — 27 000 sites 
with a random initial distribution of particles and wa- 
ter molecules. Periodic boundary conditions are used to 
eliminate boundary effects. With regard to finite-size 
effects, we have found that our results are robust with 
respect to changes in the system size. Because the sys- 
tem has a large number of degrees of freedom, a repre- 
sentative sampling of the high-dimensional phase space 
is necessary to estimate thermal averages in the equilib- 
rium state. An appropriate technique which takes into 
account the effect of statistical fluctuations on the system 
is a Monte Carlo simulation. 

In a system of fixed density (i.e. in the canonical en- 
semble), every possible configuration {n^} has the statis- 
tical weight 



Wc({rii}) 



-0H et f[{n t }] 



(17) 



where the partition function Zjq of a system of N parti- 
cles is given by Eq. 10. At equilibrium the system must 
satisfy the detailed-balance condition 



w c ({rii})P n 



Wc{{n'i})Pn 



(18) 
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where P n —> n > is the transition rate from configuration 
{rii} to a new one {n^}. The relative probability to pro- 
duce configuration {n^} from the previous one {rii} thus 
becomes the ratio of the two weights, 

r = w c({K}) = e -f}(H o!i [{n[}]-H cit [{n % }\) Qg-j 

w c {{rii}) 

and depends only on the difference in free energy between 
them. This transition probability is used in the Metropo- 
lis algorithm to generate configurations from previous 
ones. Specifically, the new configuration {n^} is accepted 
if r > 1, or if r < 1 but larger than a random number 
uniformly distributed in the interval [0, 1]. 

If the new configuration is completely different from 
the previous one, the acceptance probability is rather 
low. The method therefore sweeps randomly through the 
system considering configurations which differ from the 
previous one only by single-site exchanges of particles and 
water. After a number of thermalization sweeps, during 
which the system relaxes towards equilibrium and no ob- 
servables are calculated, the system is taken to be in equi- 
librium with only thermodynamic fluctuations present. 
Thermodynamic quantities are estimated by averaging 
over the configurations which are kept during a subse- 
quent number of steps which is sufficiently large that 
a considerable portion of the total phase space is sam- 
pled. The decorrelation time of successive configurations 
in equilibrium is found to be lower than 10 Monte Carlo 
steps (one Monte Carlo step corresponds to the consid- 
eration of every site in the system once), both in the 
coexistence phase and in the homogeneous phase. Mea- 
surements are taken only every 50 Monte Carlo steps over 
a period of 500 000 steps after 100 000 initial relaxation 
steps. This process is repeated for 10 different random 
initial configurations and the observations are averaged 
over these independent simulations. In the crystal phase 
(below) the decorrelation time of consecutive configura- 
tions is longer because of the very low temperature, and 
measurements are averaged over a larger number of in- 
dependent simulations. 

In order to find the equilibrium density for a fixed tem- 
perature and chemical potential, a grand canonical sam- 
pling (i.e. the number of particles is not constant) of the 
phase space is performed. The procedure is the same as 
in the canonical case, except that the weight of a config- 
uration {rii} in the grand canonical ensemble is 

w gc ({n t }) = , (20) 

where the grand canonical partition function 3 is given 
by Eq. (JSJ. This leads to a relative transition probability 
for configuration {n^} from a previous one {rii} which 
depends on the difference in free energy of the two con- 
figurations 

r = e -/3(fff«[{na]-#S[{»i}]). (21) 



The behavior of a system which is heated at constant 
density can be analyzed by first determining the equilib- 
rium density using grand canonical sampling at a chosen 
starting temperature and then continuously raising the 
temperature while applying a canonical sampling proce- 
dure in which the number of particles remains constant 
but particles and water may exchange sites. 

Starting from a low temperature Tq at fixed chemical 
potential p,, we observe the thermal equilibrium state of 
the system at progressively higher temperatures: for the 
same parameters as in Sees. ^ and lllll a jump in the 
density from p c \ to p c i occurs for values of p between 
p* L = 1.8 (at the lower critical temperature T£ = 0.54) 
and = 5.95 (at the upper critical temperature T v — 
2.07) (see Fig.gJ. 




P* P 
K p K p 

FIG. 5: p-T phase diagram for an aqueous solution of hy- 
drophobic particles obtained by Monte Carlo simulations, il- 
lustrating the coexistence curve. p v indicates the particle 
density. The system size is N = 27 000. The dotted line 
represents a system of density p p = 0.5 which is heated from 
a starting temperature To. At Tl a phase transition occurs 
from the homogeneous state to an aggregated phase. Further 
heating results in a second transition at Tu, where the system 
disaggregates and the particles are dissolved again. 

Fig. [SI displays the p-T phase diagram obtained from 
Monte Carlo simulations. The system shows a phase 
transition from a homogeneous state to a two-phase ag- 
gregation state at a lower transition temperature Tl , and 
a disaggregation at an upper transition temperature Tjj, 
i.e. a closed- loop coexistence curve is found. 

After rapid cooling at constant density po from a tem- 
perature Tq > Tu to a fixed temperature T < Tjj in 
the coexistence region, the system develops a clear phase 
separation into a phase with density p C 2 and nearly pure 
water. All of the particles aggregate to clusters, and after 
a certain period, one single cluster of density p c i remains, 
which occupies a fraction V C 2 = ^"'^^ of the total vol- 
ume. 

An analytical solution is possible for the ground state 
of the system at vanishing temperature, which provides 
a test for the Monte Carlo simulations. The calculation, 
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which uses the fact that E ds > and E os < 0, results in 
S(T^O) = ]im(-~ln[q as e-P E °>+q ds e-e E ««}) 

= hm(-i[ln(g os )+ln( e -^)]) 

p^oo p 

= E os (22) 

and, analogously, B(T — *• 0) = E ob . At T — 0, minimiz- 
ing the free energy is equivalent to minimizing the Hamil- 
tonian in Eq. 0. This leads to three different phases 
depending on the chemical potential fi. For fi > 2.0, the 
free energy is minimized by a solid phase (p p = 1 and 
H = 0), while for /i < —5.0 no particles are dissolved 
and pure water is found (p p — and H = B + fi). In 
between, the system forms a dispersed crystal structure 
in which the particles are arranged in such a way that 
every water molecule is the neighbor of exactly one par- 
ticle, which leads to a particle density p p — = 1/7 
and if = z{S + (x)/{z + l). 

Monte Carlo simulations of the above system at suffi- 
ciently low temperatures confirm this behavior, as shown 
in the phase diagram in Fig.0] The dot-dashed line shows 
the evolution of a system which is heated at constant par- 
ticle density p p = -^tj, which corresponds to the critical 
particle density p* found by the mean-field calculation, 
starting in the dispersed crystal phase. The critical parti- 
cle densities p* determined by mean-field approximation 
and by Monte Carlo differ slightly, but the discrepancy 
lies within the error of the Monte Carlo simulations. 



V. CHAOTROPIC EFFECT 

The addition of urea to an aqueous solution of hy- 
drophobic molecules may affect the properties of the lat- 
ter in a way that destabilizes aggregation of the non- 
polar solute In the case of protein solutions, this 
destabilization can result in a complete denaturation of 
the proteins, and even prevent their aggregation if the 
latter is due to hydrophobic interactions 33]. A high- 
concentration solution of urea is therefore often used as 
a protein denaturant. 

The underlying cause of this process, known as the 
chaotropic effect, is generally believed to be a decrease 
in the order of the water structure ('chao-trope' = dis- 
order maker), thus indirectly increasing the solubility 
of non-polar solutes [34). However, different attempts 
to discover how chaotropic agents perform this function 
have not yet been able to explain in a satisfactory man- 
ner the exact mechanism for disruption of the hydrogen 
bonds which stabilize the aggregate. In the remainder of 
this section, we adapt the MLG framework by including 
chaotropic cosolvent molecules with the aim of demon- 
strating the chaotropic effect in our model system. 

Chaotropic substances are in general those which are 
less strongly polar than water. In aqueous solutions of 
non-polar species they act to reduce the number of possi- 
ble intact hydrogen bonds between water molecules, both 



in the solvation shell and in the bulk, compared to a pure 
water-solute mixture. Within the adapted MLG frame- 
work, a straightforward approximation to the effect of 
a chaotropic cosolvent is to consider that its addition 
to strongly hydrogen-bonded, "ordered" clusters creates 
"disordered" clusters with additional broken hydrogen 
bonds and higher net energy. The creation of disordered 
states from ordered ones in the presence of a chaotropic 
cosolvent increases the degeneracy of the former at the 
expense of the latter. 

Because the coarse-grained model treats each water 
site as containing a number of water molecules, the cosol- 
vent molecules, which are generally rather small in com- 
parison with the hydrophobic particles, may be included 
implicitly at each water site by adapting only the de- 
generacies of the energy levels of the water. The states 
of water clusters containing cosolvent are thus assigned 
the degeneracies q as , u = q os - rj s and q dSjU = q ds + r] s , 
and q obtU = q ob - r} h and q dbtU = q db + rj b , where u de- 
notes urea, which we adopt as an illustrative example of 
a small chaotropic cosolvent. The cosolvent is taken to 
affect only the number of hydrogen bonds formed, and 
not their strength, so in the bimodal approximation the 
energies of the states remain unchanged. 

The effect of a chaotropic cosolvent is much stronger in 
the bulk than in the shell. While in ordered, bulk water 
both tangentially and radially oriented hydrogen bonds 
may be broken, in shell water there exist no radially ori- 
ented bonds, and the ordered bulk configurations with 
radially oriented water have already become disordered 
configurations on substitution of the central water by a 
non-polar solute particle. The number of configurations 
available to be transformed from ordered to disordered 
by the addition of cosolvent is therefore much higher in 
the bulk than in the shell. Indeed, the number of config- 
urations of ordered bulk states may be reduced almost to 
that of ordered shell water in the presence of a strongly 
chaotropic cosolvent, because the primary difference be- 
tween the two is the absence of radial hydrogen bonds 
to the central molecule in the shell of a non-polar solute 
particle. For the calculations to follow we have chosen 
the values rj b — 9.0 and rj s = 0.1, which are found to be 
suitably representative of a water/solute/cosolvent sys- 
tem. 

Within this approach the addition of a concentration 
c of cosolvent to the water leads to an additional term in 
the partition function 



ZU = £{»,} 11(9" 



-0E O 



q ds e 



x(q s,u& ^' 
x (q be~ pE ° b + q d b£ 



-AEd.Wl-AiXl-c) 



q ds ,v.e r '~" 3 Y 
pE db \ ni \i(i-c) 



Here the concentration of urea in the bulk is assumed 
to be identical to the concentration in the first sol- 
vation shell of a hydrophobic particle. On defin- 
ing B u = -ilnfeo^e^^ + q db<u e-^»] and S u = 
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— ^ ln[g OSi „e ' 3£ ' os + q<i s ,u& P Eds ], the effective Hamilto- 
nian function becomes 

N 

H&[{m}] = + J^«iAi], (24) 

i=i 

where the effective interaction J" ff and the effective chem- 
ical potential ^" ff are given by 

J c u ff = S-«S + c(B„-B + «S-.S„) (25) 
fa = S + fi + c(S u -S). (26) 

Because we are concerned with the hydrophobic interac- 
tion between solute particles, we again express the Hamil- 
tonian in terms of particles = 1 — m, 

N 

#eT[«] = K U ~ 5>>* + JeAm ~ l)Ai], (27) 

i=l 

where X" = N/i e g and A, = ~ m j)- The effec- 

tive interaction between the hydrophobic particles thus 
depends on the concentration of chaotropic agent in the 
solution. Because S u > S and B u > £>, is larger than 
yiteff and therefore it is easier to bring the non-polar solute 
into solution in the presence of chaotropic substances. In 
addition, the relation (B u — B) > (S u — S) results in an 
increase in J" ff when raising the concentration c. Hence, 
the chaotropic particles support the repulsive force be- 
tween non-polar solute molecules. The coexistence region 
of the phase diagram consequently shrinks with increas- 
ing c, and for extremely high concentrations may even 
disappear entirely. 




'0 0.2 0.4 0.6 0.8 1 
p*=l/(z+l) P P 



FIG. 6: p-T phase diagram for a ternary system consist- 
ing of water, hydrophobic particles, and chaotropic cosol- 
vent, obtained by mean-field calculation, showing coexistence 
curves for different cosolvent concentrations. The • represents 
a multi-critical point at the critical cosolvent concentration 
c* = 0.935 where the coexistence curve has shrunk to a single 
point. 

Fig. shows the effect on the coexistence curve for 
different urea concentrations c using r/b — 9 and rj s =0.1. 



The expected increase in solubility is confirmed, and the 
LCST and UCST approach each other with increasing 
cosolvent concentration. At a critical concentration c* = 
0.93, the closed-loop curve shrinks to a single point (T* — 
1.3, fi* = 2.84 and p* = l/(z + 1)), which represents a 
double critical point |35j . For concentrations higher than 
c* , the aggregation phase disappears completely. 



VI. DISCUSSION 

We have studied the aggregation of hydrophobic so- 
lutes in water using the bimodal MLG model, which we 
have adapted to describe a coarse-grained system where 
each site can be occupied by one or more molecules. One 
of the objectives was to reproduce the thermodynamic 
properties associated with the hydration and aggregation 
of non-polar solutes in water. For this purpose, Monte 
Carlo simulations were conducted to establish the phase 
diagram for the process, and the results were compared 
to a mean-held calculation. 

As expected, both methods display clear phase transi- 
tions at a LCST and an UCST within a range of densities, 
and thus define a closed-loop coexistence region. Outside 
this region the system appears as a homogeneous particle- 
solvent mixture, while inside it a separation occurs into 
two phases of fixed (upper and lower) coexistence den- 
sities. The exact time evolution of the system after a 
quench into both the metastable and the spinodal regions 
at constant density has not yet been studied, but would 
offer interes ting insight into liquid-liquid phase-ordering 
kinetics [Mill. 

At low temperatures (To in Fig. |SJ) the system is in 
the homogeneous region where the solubility of the so- 
lute is high and the solvent-induced effective interaction 
between solute particles repulsive. On raising the tem- 
perature at constant density, the system shows a sharp 
transition to an aggregation state at the LCST Tl. In the 
equilibrium state at this temperature, clusters of aggre- 
gated, hydrophobic solute molecules with density p c2 are 
suspended in nearly pure water (density p c \ ~ 0). The 
aggregation allows the solute particles to minimize their 
exposed surface and to reduce the structural enhance- 
ment of the surrounding water, thus causing a positive 
entropy change. The effective interaction between non- 
polar particles becomes attractive above Tl, and there- 
fore the formation of aggregates is preferred over the so- 
lution of single solute particles. Further heating of the 
system results in a second phase transition at the temper- 
ature T\j where the hydrophobic particles disaggregate, 
and above this temperature one homogeneous phase is 
found. This process is dominated by the favorable en- 
tropy change of solvation in the whole system. 

As can be seen from the \i-T phase diagram in Fig. 0] 
the transition line for the Monte Carlo simulations and 
the mean-field calculation correspond rather well. They 
agree closely on the lower critical point, as well as on 
the densities above and below the transition tempera- 
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tures. The upper critical point determined by Monte 
Carlo, however, lies at a lower temperature than the 
mean-field one. This result is not surprising, because 
it is well known that mean-field calculations, which ne- 
glect fluctuation effects, generally overestimate transition 
temperatures. The good agreement at low temperatures 
is rather a signature of the predominance of local effects 
which do not involve large fluctuations from site to site. 

Apart from this difference, the coexistence curves from 
Monte Carlo and the mean-field calculations are qualita- 
tively similar. At very low temperatures, however, the 
former show an additional crystal phase which cannot 
be explained by the latter, because it arises from spa- 
tial ordering of the hydrophobic solute and is therefore 
neglected by the present mean-field considerations, al- 
though this phase could be recovered within a more re- 
fined mean-field approximation. The appearance of this 
phase confirms the prediction of our analytical calcula- 
tion at T = 0. 



250 
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FIG. 7: Closed-loop solubility curve for nicotine/water, cal- 
culated within the mean-field framework by using the experi- 
mental value of the critical density and renormalizing by com- 
parison with Monte Carlo simulations (see text). The curve 
shows good quantitative agreement with the experimentally 
determined solubility curve (o) reproduced from Ref. 0. 

Encouraged by the qualitative and quantitative agree- 
ment of the mean- field calculation and Monte Carlo simu- 
lations, we have taken experimental values for the critical 
density p* in different systems to adapt the parameters 
of our model. Analytically, we have obtained the relation 
p* = j-jtj by mean-field calculations, and have confirmed 
this numerically by Monte Carlo simulations. Thus from 
an experimental value for p* we may extract an effective 
coordination number z cS , which is then introduced in the 
calculations. The critical density for the system nico- 
tine/water is p* = 0.4 ;7J, which results in z^ w = 1.5. 

This coordination number z cS may be interpreted as 
the average number of hydrophobic solute molecules sur- 
rounding any chosen solute particle, which is relevant for 
the net effective hydrophobic interaction leading to at- 



traction or repulsion. After changing the values for the 
energy levels to E ds = 3.4, E db = 3.3, E ob = —3.3, and 
E os = —4.1, the closed-loop coexistence curve is in good 
agreement with the experimental curve (Fig. [7J. The 
UCST is higher than the measured one, as expected be- 
cause of the mean-field nature of the calculation. The 
ratio between the mean-field result Tu^i and the exper- 
imental value Tj/ jeX p is Tu, m {/Tu,cxp = 1-25. In fact this 
value agrees quantitatively with the ratio Tu, m {/Tu,MC = 
1.35 which we obtain by comparing the mean-field calcu- 
lation with the more accurate Monte Carlo simulations 
when both are performed for z — 6 (cf. Fig. 0J. Us- 
ing this as an effective scaling factor to renormalize the 
mean-field results for different z cS yields good agreement 
with the experimental results for the nicotine/water sys- 
tem (Fig. UJ . We have repeated this procedure for the 
system poly(ethylene glycol) in water, which has a crit- 
ical density of p* = 0.15 and a much larger molecular 
weight M w = 3350 8], and which resulted in a similar 
agreement with the experimental curve. 

In addition to studying hydrophobic solutes in pure 
water, we have analyzed the effect of a chaotropic sub- 
stance in the solution. The destabilizing effect can be at- 
tributed indirectly to a disordering of water molecules by 
chaotropic cosolvents associated with a rearrangement of 
the energy states of the water. As a first approximation, 
we have included the chaotropic agent in liquid water by 
increasing the number of disordered water states with re- 
spect to ordered ones. A mean-field calculation for the 
model with cosolvent shows the expected increase in sol- 
ubility of the hydrophobic particles. At concentrations 
higher than a critical value, the aggregation phase disap- 
pears completely. Monte Carlo simulations confirm this 
behavior. 

Experimentally, high concentrations of chaotropic sub- 
stances are required to cause destabilization of aggregates 
of hydrophobic particles. As an example, the solubil- 
ity of the highly hydrophobic amino-acid phenylalanin at 
room temperature is doubled in an 8-molar urea solution, 
which consists of an equivalent volume fraction of urea 
and water. While within the current framework it is diffi- 
cult to express the cosolvent concentration explicitly be- 
cause the model contains more than one water molecule 
per site, the results are qualitatively correct. A cosol- 
vent concentration of at least 50% is required to reduce 
significantly the extent of the solubility region (Fig.EJ. 

We should emphasize that in the mean-field analysis 
presented here the concentrations of chaotropic agents in 
the bulk and in the first solvation shell were assumed to 
be equal. However, chaotropic molecules are found pref- 
erentially in the solvation shell of non-polar solute par- 
ticles. Considering two different concentrations may be 
expected to provide better insight into the exact mecha- 
nism of the chaotropic effect |2g . 

Overall, we have shown that qualitative features of the 
liquid-liquid demixing process of hydrophobic aggrega- 
tion, as well as of the chaotropic effect, may be explained 
satisfactorily within a simple model for aqueous solutions 
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of non-polar particles by including hydrophobic interac- 
tions only in terms of changes in water structure. Al- 
though the explicit terms of model describe solely the 
states of water molecules in solution, we have demon- 
strated that it contains implicitly both two- and even 
many-particle interactions between hydrophobic solute 
molecules. The complete density-temperature phase dia- 
gram was established, by both analytical and numerical 
techniques, and illustrates the characteristic properties 
of hydrophobic aggregation. 
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